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The analysis of polycrystalline materials benefits greatly from accurate quantitative 
descriptions of their grain structures. Laguerre tessellations approximate such grain 
structures very well. However, it is a quite challenging problem to fit a Laguerre tes¬ 
sellation to tomographic data, as a high-dimensional optimization problem with many 
local minima must be solved. In this paper, we formulate a version of this optimization 
problem that can be solved quickly using the cross-entropy method, a robust stochastic 
optimization technique that can avoid becoming trapped in local minima. We demon¬ 
strate the effectiveness of our approach by applying it to both artificially generated and 
experimentally produced tomographic data. 
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1. Introduction 

In recent years there have been significant advances in the tomographic characteri¬ 
zation of materials. As a result, it is now possible to carry out detailed investigations 
of the 3D grain structures of polycrystalline materials; see, e.g., m \ • A fundamen¬ 
tal ingredient in any such investigation is a suitable quantitative description of the 
grain morphology. Such a description contains the key features of the structure, 
ideally free from noise and imaging artifacts. A good description usually results in 
significant data compression, describing large 3D voxel data sets using only a small 
number of parameters. Data compression is necessary, for example, when carrying 
out analysis of sequences of tomographic data sets (e.g., the high time resolution 
in-situ synchrotron images considered in 0) In addition, the description of tomo¬ 
graphic data via tessellations provides a basis for the statistical analysis of grain 
structures and, in some cases, can be used to develop stochastic models of material 
microstructures; see, e.g., IMS]- 

The most commonly used quantitative descriptions of space-filling grain ensembles 
are based on tessellations, which divide the space into disjoint regions called cells. 
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The cells represent the individual grains. The most widely used tessellation model 
is the Voronoi tessellation (see, e.g., nans!), which takes, as parameters, a set 
of generating points. The space is then divided into convex cells by assigning each 
point to its nearest generator. The Laguerre tessellation (see, e.g., [Hi US]) is a 
generalization of the Voronoi tessellation that also partitions the space into convex 
cells. Like the Voronoi tessellation, the Laguerre tessellation is generated by a set of 
points; however, unlike the Voronoi tessellation, these points are weighted, with the 
weights influencing the size of the cells. Consequently, the Laguerre tessellation is 
able to describe a wider range of structures than the Voronoi tessellation. For this 
reason, the Laguerre tessellation is a popular choice for modeling polycrystalline 
grain structures [8E1H2H22! and other materials, such as foams mmm- 

In order to describe a tessellation by a set of generating points, it is necessary 
to solve an inverse problem: that is, a set of generating points that produce the 
observed cells must be found. The Voronoi inverse problem (VIP) is well-studied; 
see, for example, [251 - 130] . Recently, Duan et al. [31] proposed an algorithm that finds 
solutions to the Laguerre inverse problem (LIP). Although the examples considered 
in [3T] are restricted to 2D, the methodology is easily applied in higher dimensions. 

The solutions to the VIP and the LIP assume that the empirical data constitute 
perfect descriptions of the observed cells. However, this is not true when working 
with tomographic data, which is distorted by noise and also contains imprecision 
arising from discretization during the imaging process. It is also worth noting that 
real-world materials are not perfectly described by Laguerre tessellations (though 
the descriptions can be quite good). These limitations mean that methods that 
attempt to invert a tessellation extracted from the tomographic data do not, in 
general, result in good fits. The LIP is solved by iteratively finding the generating 
points of the given tessellations. When applied to imperfect data, this iterative pro¬ 
cedure propagates errors, resulting in tessellations that have little correspondence to 
the tomographic data. Thus, when dealing with empirical data, it is not appropriate 
to attempt to solve the LIP. Instead, the generating points of a Laguerre tessellation 
that is a good approximation of the material must be found. This is, at its core, 
an optimization problem. We call this problem the Laguerre approximation problem 
(LAP). The corresponding Voronoi approximation problem has been considered in 
the literature, beginning with [32] . 

A simple heuristic approach for solving the LAP was proposed in [6|. More so¬ 
phisticated approaches, which formulate and solve an optimization problem, are de¬ 
scribed in [331 - 135] . Although these techniques provide good fits in certain settings, 
they are either limited to small sample sizes or require the considered tessellations 
to be sufficiently regular. 

In this paper, we present a fast and robust algorithm for fitting Laguerre approx¬ 
imations to large data sets. More precisely, we formulate an optimization problem 
where we minimize the discrepancy between the grain boundaries observed in the 
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image data and the grain boundaries produced by our Laguerre approximation. 
The cost function is chosen so that it can be evaluated very efficiently and that all 
necessary information can be easily obtained from image data. We then solve the 
optimization problem using the cross-entropy (CE) method [361139] . a stochastic op¬ 
timization algorithm that is able to escape local minima. We carry out experiments 
on both real and artificially-generated image data that show our approach is able 
to produce very good fits. 

This paper is structured as follows. In Section [2] we review some key properties 
of Laguerre tessellations. In Section [3] we give a more complete description of the 
LAP and formulate our optimization problem. Then, in Section [4] we introduce the 
CE method as a robust tool for solving this optimization problem. Section [5] gives 
results for both artificial and experimental data that demonstrate the effectiveness 
of our approach. Finally, Section [6] summarizes our results and suggests directions 
for further research. 


2. The Laguerre tessellation 

In the following section, we define Voronoi and Laguerre tessellations and give a 
number of properties that we will use to solve the LAP. For notational convenience, 
we only consider tessellations in M 3 . However, our methods are easily applied in 
other dimensions. 

The Voronoi tessellation is defined by a locally finite set of generating points, 
{xj}j g x, where ICN denotes the index set consisting of natural numbers. The cell 
corresponding to the ith generating point, x* e M 3 , is given by 

Ci = {y G M 3 : ||y - x*]| < ||y - xy || for all j G 1} , (1) 

where || • || is the Euclidean norm on M 3 . In the case of the Laguerre tessellation, 
the generating points, {(xj, rj)}j £ x, are weighted, where we assume that r* > 0. The 
cells of the Laguerre tessellation, {Ci\i & x, are then defined by 

Ci = {y G M 3 : pow(y, (x*, r*)) < pow(y, ( Xj,rj )) for all j £l}, 

for all i £ I, where the Euclidean norm used in 0 is replaced by the so-called 
power distance 


pow(y, (x,r)) = ||y - x|| 2 - r 2 . 

Given the generator points, the faces of the cells can be computed efficiently; see, 
e.g., [401 - 142] . As we take the weights r,; to be positive real numbers, the generat¬ 
ing points have a geometric interpretation as spheres. That is, we can represent a 
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generating point, (x,,rj), as a sphere with center Xj and radius r*. 

The flexibility of the Laguerre tessellation comes at a cost. For example, while 
each cell of a Voronoi tessellation contains its generating point, the generating points 
of a Laguerre tessellation may not be contained in their corresponding cells. In some 
cases, a generating point of a Laguerre tessellation may not even produce a cell. That 
is, it is possible that there exists a point, (xj, r*), such that C* = 0; see, for example, 
m- In addition, while the cells of a Voronoi tessellation uniquely determine its 
generating points, there are uncountably many sets of generating points that can 
generate a given Laguerre tessellation; see, [31, Mi S3]- Thus, while the VIP has a 
unique solution, the LIP has uncountably many solutions. Under mild conditions, 
however, it can be shown that the generating points of a Laguerre tessellation are 
uniquely determined given one generating point and the weight (or one coordinate) 
of a generating point in an adjacent cell; see m- 

We will often find it convenient to consider planes that are equidistant from 
two generating points (under their respective power distances). That is, given two 
generating points, (x,;, r t ) and (xj , r 2 ), we consider the separating plane given by 

Pi,j = {y G M 3 : nj d y + b i;j = 0} , 


where 


n i,j 


2 • (xj - Xj) 

2 ■ ( Xj -Xi)|| 


is the unit normal vector of Pij and 


||x t || 2 - 11 Xy 11 2 + r 2 j - rf 
1,3 || 2 - (xj -Xi)|| 

is the shortest distance from the origin to the plane; see, e.g., [43; Section 2.1] for 
more details. The plane Pij defines a half-space 

Hij = {y £ K 3 : njj-y + b itj < 0} , 

which covers the cell C%. Note that the intersection of all such half-spaces, 
f1fce2\{i} defines the cell C\. Because the normal vectors are taken to be unit 
vectors, the distance from an arbitrary point, x G M 3 , to the plane Pij is given by 

d(x,Pij) = |nJjX + &ij|. (2) 
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3. The Laguerre approximation problem (LAP) 

Given an exact description of a Laguerre tessellation (e.g., in terms of the half-spaces 
defined above), it is not too difficult to solve the LIP. That is, it is straightforward 
to find a set of weighted points that is able to generate the given tessellation. 
Furthermore, these points can be chosen to satisfy certain constraints; see ED- 
When dealing with tomographic data, however, the description of the tessellation is 
not exact. This is because, even if the material itself can be perfectly described by 
a Laguerre tessellation, the noise and discretization errors inherent in the imaging 
process mean that the cell boundaries extracted from the data are subject to error. 
The LIP can be solved as described in ED- However, the iterative computation of 
generators is sensitive to imperfect data and errors propagate quickly. As a result, 
the ensuing tessellations often do not correspond well to the tomographic data. 
Therefore, solving the LIP for empirical data is generally ill advised. Instead, we 
wish to find generating points that produce a Laguerre tessellation that is as close as 
possible to the tomographic data (with respect to some metric). This is, at its core, 
an optimization problem. In order to properly formulate this optimization problem, 
a discrepancy measure needs to be defined. We then choose the generating points 
of the approximating tessellation in order to minimize this discrepancy. The choice 
of discrepancy measure depends on the nature of empirical data. In this paper, we 
work with tomographic data. 

3.1. Tomographic image data 

We assume that the tomographic image data constitutes a collection of voxels in a 
convex window. Furthermore, we assume that the image has already been segmented 
— e.g., by the watershed transform; see }451 06j — and that the voxels have been 
labeled by their corresponding grains. Thus, the data is of the form {/(x,y, z) £ 
{0, ..., N } : (x, y, z) E W}, where N > 1 denotes the number of grains and W C N 3 
is a grid of voxel coordinates. The ith grain region, Ri, is then given by Ri = 
{(x,y,z) 6 W : I{x,y,z ) = i}. Note that Ro does not correspond to an actual 
grain, but is either the empty set or a collection of one voxel thick layers that 
separate grains. Such thin layers often arise in segmentation procedures such as the 
watershed transformation. A grain region Ri may correspond to the grain itself or 
a region that contains the grain (if the space is not completely filled with grains), 
cf. Figure [lj We assume that the grain regions are roughly convex and that the 
segmentation is of a high quality. 

3.2. Discrepancy measures 

The most direct way to measure the discrepancy between 3D image data and a 
Laguerre tessellation generated by the points {(x,;, rj)}j £ x is by counting the number 
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Figure 1. Illustration of voxels labeled by their corresponding grains. The grain regions are assumed to form 
a tessellation with convex cells. The grains themselves may not fill their corresponding regions completely, as 
depicted for grain number 3. Here, we show the grain region (i.e., voxels with label 3), a cell approximating 
the grain region (colored cyan) and the grain itself (colored blue). 


of incorrectly assigned voxels. That is, we calculate 

^vox ({(xiXhex) = # | (x,y,z) <E W : I(x,y,z ) > 0 ,I(x,y,z) A I(x,y,z) j , 

where I is a discretized version of the tessellation generated by {(xj,rj)}iex and 
A A denotes the number of elements in some set A. Many other discrepancy mea¬ 
sures considered in the literature are of a similar form. For example, in [35], the 
difference between a cell in the approximating tessellation and its equivalent in the 
empirical data is measured by the intersection of the approximating cell with the 
corresponding adjacent cells in the empirical data. Minimizing the total overlap of 
the approximating cells is equivalent to minimizing the number of incorrectly as¬ 
signed voxels. We call discrepancy measures that aim to minimize the volume of 
the difference between the empirical and approximating tessellations volume-based, 
measures. 

Although minimizing a volume-based measure gives a very good fit to the data, 
such measures are poorly suited to most numerical optimization methods. This is 
because evaluating such a discrepancy measure is computationally expensive. For 
example, when considering the number of incorrectly assigned voxels, a discretized 
version of the approximating tessellation needs to be generated. In practice, this 
limits the number of times such a discrepancy measure can be evaluated. As a 
result, approaches that aim to minimize volume-based discrepancies are restricted 
to small data sets with small numbers of grains (e.g., 109 grains in |34| ) or are forced 
to use non-optimal optimization techniques, such as gradient-descent, which do not 
require too many evaluations of the discrepancy and, as such, limit the time taken 
by the fitting procedure (ideally to substantially less than 24 hours). 

As we propose to use stochastic optimization techniques to solve this problem, we 
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need a discrepancy measure that is significantly faster to evaluate. We can achieve 
this by considering an interface-based discrepancy measure — a measure that con¬ 
siders only the boundaries between cells — instead of a volume-based one. If an 
approximating tessellation can accurately reproduce these interfaces, it will be a 
good approximation of the data. The primary advantage of interface-based discrep¬ 
ancy measures is that they can be calculated from the generating points of the 
approximating tessellation without the need to generate the tessellation itself. 

In order to define such a discrepancy measure, we consider the sets of voxels that 
separate adjacent cells. We define the interface between two grains, i and j, in the 
empirical data by 


N it j = {(x,y,z) eW : Af 2 e(x, y,z)r\Ri^tt and Af 26 (x, y, z ) n Rj / 0} , 


where Af 2 e(x, y, z) denotes the 26-neighborhood of ( x,y,z ) — i.e., the voxels 
(x',y',z') that have a distance less than or equal to \/3 from ( x,y,z ). This set 
contains all voxels that touch both grains. Note that Nij = 0 if the grains are not 
adjacent. If the tessellation generated by {(xj,ri)}j e x is a good approximation to 
the empirical data, then the plane separating the generating points of adjacent cells 
i and j (which defines the boundary of the two cells) should be close to IVy. Thus, 
we measure the distance between Wj and P t .j, the plane separating cells i and j in 
the approximating tessellation. The discrepancy between the approximating tessel¬ 
lation and the empirical data is then given by the sum of squares of these distances. 
That is, we define 


N N 

2? ({(xi,r i )}iex) = 5^ 2 (3) 

j= i k=j+ 1 yeiV,, fc 

where d( y, Pj.k) is given in ([2]) and {Pj.k}j,kei is the set of separating planes de¬ 
termined by the generating points {(x.j, ri)}i £ x- Note that this discrepancy measure 
can be calculated without generating the approximating tessellation. A separating 
plane can be computed for all pairs of generators (or generator candidates, when 
performing the optimization) — even if their cells are empty or not adjacent in 
the Laguerre tessellation. Therefore, such ‘degenerate’ cases are not a problem for 
our approach. By matching all separating planes to the corresponding test points 
at once, Laguerre cells computed based on ‘good’ configurations of generators will 
match the desired cell configuration well enough. 
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3 . 3 . Minimizing the discrepancy 


When aiming to minimize a volume-based discrepancy, the LAP reduces to the 
optimization problem of finding 


{( x ;,U*)}iex = argmin P vox ({(xi,r-i)} ie x). (4) 

There are a number of significant difficulties which must be overcome in solving Q . 
In particular, 

(i) the optimization problem is high-dimensional, 

(ii) the optimization problem has many local minima, 

(iii) the discrepancy is expensive to evaluate. 

The approach developed in [33] avoids problems (ii) and (iii) by finding the exact 
solution of a linear program (an optimization problem where a linear cost function 
is minimized subject to linear constraints). However, when the LAP is transformed 
into a linear program, the size of the resulting problem grows very quickly in both the 
number of voxels and the number of grains. This means that the linear programming 
approach cannot be applied to large 3D data sets containing many grains. 

In [35j, gradient-descent methods are used to obtain Laguerre approximations 
(with a relatively low number of evaluations of the discrepancy). Unlike the linear 
programming approach, the dimensionality of the optimization problem grows only 
linearly in the number of grains (and does not depend on the number of voxels). 
Thus, [35] avoids problem (iii) and, to a lesser extent, problem (i). However, the 
gradient-descent methods used there become stuck in local minima. Thus, the ap¬ 
proach considered in [35] is reliant on the ability to find initial conditions that are 
close to global optima. Although it is often possible to find good initial conditions 

- e.g., when fitting tessellations to foams with regular structures — this is not 
always the case, as will be seen in Section [5] In addition, even when good initial 
conditions can be found, the quality of the approximating tessellation can usually 
be significantly improved when the optimization algorithm is able to escape local 
minima. 

Stochastic optimization methods are widely used tools for solving high- 
dimensional optimization problems with many local minima; see [38] for an overview. 
Such methods have been used to solve problems related to the LAP. In [31], the 
CE method is used to find a solution to the LIP and, in [33], simulated annealing 
is used to fit 3D Laguerre tessellations to 2D image data. However, no stochastic 
approach has yet been proposed in order to solve the LAP for large, noisy 3D data 
sets. 

In our approach to solving the LAP for large tomographic data sets, we use the CE 
method to minimize an interface-based discrepancy measure. The CE method has 
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been widely applied to high-dimensional multi-extremal optimization problems; see, 
(36b39j . It has a number of advantages over other stochastic optimization methods. 
For example, simulated annealing (described in [38102]) cannot be easily applied to 
the LAP, as it is very difficult to find an appropriate cooling schedule. In addition, 
the CE method can be easily parallelized. 

By using the CE method to minimize an interface-based discrepancy measure, 
we have a method that can escape local minima and solves a problem whose di¬ 
mensionality grows linearly in the number of grains. Although the discrepancy is 
expensive to evaluate, we reduce this cost significantly by minimizing an interface- 
based measure. In addition, we further reduce the cost of calculating the discrepancy 
by replacing T> defined in ^ with an approximation, T>. 


3.4. Approximating the discrepancy 

The interface-based discrepancy, T>, measures the distance between voxels on the 
boundaries between cells in the empirical data and the planes separating the gener¬ 
ating points in the approximating tessellation. In order to calculate this discrepancy, 
every voxel in each grain boundary is considered. Although this is much faster to 
calculate than a volume-based discrepancy, it is still computationally expensive. 
However, very good approximations of T> can be obtained by instead considering 
sets of test points that describe the interfaces between the grains. We find that 
test points obtained by fitting approximating planes to the interfaces, using the 
orthogonal regression approach introduced in jT8j, work very well. 

In order to calculate the test points, we consider only the points separating the 
two cells being considered (i.e., we ignore points of contact between three or more 
grains). Thus, instead of Nij. we consider 

N tj = {( x iV,z) € N i.j ■ ( x ,V,z) i N i,k and ( x,y,z ) ^ N kJ ,k € {l,...,N},k £ {i,j}}. 

This avoids some numerical issues that could arise when using orthogonal regression. 

We then approximate the boundary between the zth and jth grain using the plane 
that minimizes the total squared distance to all points in I\W. Determining this 
plane is a least-squares problem, which we solve using singular value decomposition. 
The approximating plane passes through the centroid, c 6 M 3 , of the set N*y 
We obtain a normal vector, n E M 3 , to the plane by taking the right-singular 
vector corresponding to the smallest singular value of the matrix containing the 
voxel coordinates in IV*- shifted by —c. For more details on the singular value 
decomposition approach, see [49] . The test points, C M 3 are chosen to lie on 
this approximating plane. More precisely, we put a circle with radius yj'ftN *■/4 
around the centroid c with the same orientation as the plane. We then place 10 test 
points equidistantly on this circle. An illustration is given in Figure [2] 
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circle located on fitted plane 



Figure 2. Illustration of test points construction. First, a circle is fitted to a grain interface. The center 
of the circle is given by the centroid of the voxels adjacent to both grains, its 3D orientation is obtained 
from the orthogonal regression plane. The radius is chosen with respect to the interface area. Finally, 10 
test points are placed equidistantly on this circle. 


In cases where the estimation of the plane is not ‘stable’ (i.e., the number of 
voxels in IV*- is too low to determine the location and orientation), we simply use 
the centroid c for all 10 test points. The criterion for ‘stable’ is simple: the smallest 
singular value must be smaller than half the second largest one. This way we can 
be confident that the third and shortest axis (which is perpendicular to the plane) 
is correctly identified. 

Given the test points for the boundary between adjacent grains i and j, we define 
the approximate discrepancy at the boundary by 

Vij (xj, n, Xj, rj) = d( y, P itj ) 2 , 

ysT.j 

where d(y, Pj,k) is given in ([ 2 ]) and P ltJ is the plane separating the generating points 
(xj,7-j) and (xj, rj ). The total approximate discrepancy is then given by 

1 N N 

P({(xi,ri)} ie i) = — w N -— J2 ^j,k(^j, r j,^k,r k ). (5) 

2—tj=i 2^k=j+i tF j,k j = i k=j- i-i 

Note that, here, we normalize the discrepancy. This is done in order to make the 
cost function easier to interpret: it can be thought of as the average squared distance 
of test points from their corresponding separating plane. 


4. Solving the LAP using the CE method 
4.1. The CE method 

The cross-entropy (CE) method is a stochastic optimization method that is able 
to solve many difficult optimization problems, including combinatorial optimiza¬ 
tion problems and continuous optimization problems with many local minima; see 
jSH [39] . The fundamental idea of the method is to describe the location of the 
global minimum of an m-dimensional cost function, S(-), in terms of a degener- 
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ate m-dimensional probability distribution. That is, a random variable with this 
distribution takes only a single value, namely x* = argmin xgRm S'(x). If there is 
more than one global minimum, each minimum will have a corresponding degen¬ 
erate distribution. The CE method works by ‘learning’ one of these distributions. 
In the continuous setting, this is done as follows. A parametric density, /(x;0), is 
used to describe the possible locations of a global minimum. A sample of size M is 
drawn from this density and ordered by the corresponding values of S. The \pM~\ 
(where ["•] is the ceiling function) members of the sample with the lowest values 
of S are identified as the ‘elite’ sample, where p E (0,1). The elite sample is then 
used to update the parameter vector, 0. The updating is done by choosing 6 to 
minimize the cross-entropy distance (also called the Kullback-Leibler divergence) 
between /(x; 6) and the targeted degenerate probability distribution. The process 
is continued until a stopping condition is met, e.g., the probability distribution is 
nearly degenerate. Thus, the CE method uses information about good choices of 
arguments to find better choices. The CE method and its convergence properties 
are discussed extensively in |36l - f39] . 

Usually, the parametric density, /(x;0), is chosen to be a product of normal 
densities. That is, we choose 

/(x; 0 ) = (p(xi; p!, a 1 )(p(x 2 ; P2, CT2) ■ ■ ■ l^rm ®m) •> 

where ip(-;p,a) is the density of a normal distribution with mean p and standard 
deviation a and 6 = Thus, each component of the argument 

of S(-) is drawn independently from a normal distribution. There are two main 
reasons for using normal distributions. First, as the standard deviation of a normal 
distribution goes to zero, the distribution converges to a degenerate distribution. 
The second reason is that the parameter vector, 9*, which minimizes the cross¬ 
entropy distance, is simply given by the maximum-likelihood estimates obtained 
from the elite sample. In other words, we update the parameters by setting the 
means equal to the sample means of the elite sample and the standard deviations 
equal to the sample standard deviations of the elite sample. 

Using this approach, we have the following general algorithm for estimating x* = 
arg min xgRm S(x). 

1) Initialization. Identify an initial parameter vector, 

0 (o) = (/4°\ Oi°\ • ■ •, and set n = 0. 

2) Sampling. Draw an independent sample X^ 1 ),..., X^'U from /(x;£^ n )) and 
sort it so that 5(X (1 )) < ■ ■ ■ < S^X^U). We denote the jth component of X^ 

by xf. 

3) Updating. Calculate the sample means and standard deviations of the elite 
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sample. That is, calculate 
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4) Iteration. If a predetermined stopping condition is met, terminate. Otherwise, 
set n = n + 1 and repeat from step 2. 


4-1.1. Parameter choice and stopping conditions 

The parameters of the CE algorithm help to improve the speed of convergence and 
quality of solutions. For example, a good choice of 6^ improves the convergence 
properties of the algorithm. The means, ■ ■ ■, pin , should be chosen as close 
as possible to the optimal solution. The initial standard deviations, ..., ain , 
should be chosen large enough that the algorithm is able to escape local minima 
but not so large that good configurations are quickly abandoned. 

Ideally, the size of the sample in each step, M, should be quite large. However, 
there are often memory and performance constraints that limit this size. The pa¬ 
rameter p, which controls the size of the elite sample, should be chosen large enough 
that a representative sample of good solutions are included in the elite sample. How¬ 
ever, if p is chosen too large, the algorithm will take too long to converge. In general, 
the bigger M is, the smaller p should be. 

A standard stopping condition for the algorithm is that it terminates when the 
cost function does not decrease significantly for a given number of steps. 

4-1-2. Variance injection and dynamic smoothing 

When using the CE approach outlined above, the standard deviations of the normal 
densities may shrink too quickly. In this case, the CE algorithm can converge to a 
sub-optimal solution. In order to guard against this, we use variance injection. The 
basic idea is to occasionally increase the variances of the distributions so that the 
algorithm can easily escape local minima. Usually, variance is injected when the cost 
function does not decrease significantly enough over a given number of iterations. 
The magnitude of this increase can depend on the current value of the cost function. 
If variance injection is used a number of times without a significant decrease in the 
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cost function, then the algorithm is terminated. 

An alternative to variance injection is to use smoothing when updating 6. The 
updating step for the means at the nth iteration is then of the form 

^ +1) =a-X 1 + (l-a)-^\ ..., + 

and the updating step for the standard deviations is given by 

=/?•<?! + ( l-/3)-<r[ n \ a^=P-S m + (l-P)-a^, 

where a, (3 G (0,1] (typically [0.7,1]). It is also possible to carry out dynamic 
smoothing, with both a and (3 taken to be functions of the number of steps. In 
our experience, variance injection is more effective than smoothing. For a detailed 
discussion of both variance injection and dynamic smoothing, see m- 


4.2. Solving the LAP 


We solve the LAP by minimizing the approximate discrepancy described in Section 


3.4 That is, we solve 


{( x *Vi)}*ex= argmin V ({(x*, ri)} i£X ) 

{(xi.rdhex 

N N 

= arg min E E v j,k (xj, r j,*k,rk), (6) 

{(xordbex j=1 k=j+1 


where V is given in ([5]). Putting this in the terminology of Section 
arguments that minimize the cost function D , namely the N generating points, each 
of which is described by three coordinates and an associated radius. We associate 
a normal distribution with each of the m = 41V values that need to be determined. 
Thus, for i = 1 the coordinates and radius of the zth generating point, 

(xj,rj) = Zi,ri), are each described by a normal density. We denote the 

initial means and standard deviations by fi'x }, g'y }, , h'r} and , cr$. 

In order to ensure that the radii are positive, we truncate the corresponding normal 
densities to the positive real line (this is possible without changing the updating 
rules, see, e.g., mm)- 


4.1 


we find the 


4-2.1. Variance injection and stopping conditions 

We apply variance injection when the cost function does not decrease significantly 
over a period of r > 0 iterations. More precisely, at each iteration of the algorithm, 

~(n) 

we record , the minimal value of the approximate discrepancies calculated from 
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the sample in that step. If, at the nth step, 


-T)( n ) _ 

^min 


max tg{n _ Tj n _ 1} 


V 


(n) 


<5; 


injects 


we perform variance injection. Because many generating points may already be close 
to their optimal positions and sizes, variance injection is carried out locally with a 
magnitude controlled by a parameter n > 0. This is done by calculating the local 
cost of each cell and increasing the variance of the associated densities accordingly. 
We calculate the average cost of the ith cell by 


Vi 


N 


E N JJTH Z- J 

j=llr- L i,j j=i 


T V; 


d”> 




The local cost of a cell, T>*, is defined to be the maximum of its own average cost 
and the average cost of its adjacent cells. That is, 


V* = max { T>k : k = i or k such that N* k / 0} , i = 1,..., N . (7) 


The variance injection is performed by setting 

4? ) =4? ) + «-V / W ! o$>=*«> + K. y /T>;, 

*W=a(? + K .yvf, 


for i = 1,... ,N. 

If, at any stage, the benefit of variance injection becomes negligible, we stop 

(n) 

performing it. More precisely, if the current minimum cost, PW divided by the 
minimum cost immediately prior to the last variance injection is larger than 7 E 
( 0 , 1 ), we no longer carry out variance injection. 

The algorithm is terminated when 


'T'db 

p min - maxtein—T, 2 W n 


V 


(n) 


< h 


terminate^ 


where (^terminate <C ^inject* 

4-2.2. Initial configuration 

Reasonably good initial tessellations can be obtained directly from the tomographic 
data [ 6 ]. The centroids and equivalent radii of the grains in the tomographic data 
are used for the coordinates and radii of the generating points. The centroids, 
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{(cxi,Cy.,c z J}iei, are given by 


— 


E 






C?/- — 


E 


2/ 






E 


( x,y,z)£Ri ( x,y,z)£Ri 

for z = 1,..., N. The equivalent radii are given by 


(x,y,z)eR, 


#Ri' 


n - \l 


for i = 1 ,N. The initial means of the densities are then given by = c Xi , 
fly} = c Vi , i^i = c Zi and = ri for i = 1 ,,N. 

The initial standard deviations are chosen proportional to the local cost of the 


cells in this approximation. The local costs are calculated as in Section 4.2.1 Thus, 


4? = d 0) = d“> = 4 0) = 


= sM, 


for i = 1,... ,1V. 

4-2.3. Choice of control parameters for the CE method 

Our investigations showed that a parameter choice of M = 4000, p = 0.05, 
^inject = 0.05, terminate = 0.01, T = 10, 7 = 0.9 and K = 0.25 Works Well for 
a large number of different data sets. The elite set then consists of \pM~\ = 200 
samples, which results in a sufficiently large sample of good solutions. Note that 
other parameter values may improve the speed of convergence (or, in the case of 
the parameters controlling variance injection and termination, improve the quality 
of the approximation). However, the basic performance of the CE algorithm is rel¬ 
atively robust to parameter choice. That is, the convergence behavior and quality 
of approximations are good for a wide range of parameters. 

4 . 2 . 4 . Edge cells 

The approximations obtained using the CE algorithm have unbounded cells at the 
edge of the observation window. This is because the corresponding grains in the 
tomographic data are not delimited by other grains. The unbounded cells are then 
intersected with the observation window, W. Note that, in most cases, edge grains 
are of little scientific interest and are not considered when analyzing the data. In 
practice, edge grains may be only partially observed or may not be representative 
(e.g., when stuying the dynamics of grains undergoing grain coarsening). In this 
paper, we explicitly ignore grains and their approximating cells if they lie at the 
edge of the observation window. 
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5. Experimental results 


In this section, we present the results of a number of numerical experiments which 
we carried out in order to demonstrate the effectiveness of the CE approach. We 
consider three distinct data sets. Two of these data sets are produced using stochas¬ 
tic models. These ‘artificial’ data sets allow us to investigate how well our method is 
able to reconstruct tessellations when we are certain the underlying tessellation is, 
indeed, Laguerre. The first artificial data set is produced using a stochastic model 
that describes a polycrystalline material undergoing grain coarsening. The second 
artificial data set is, by design, much more pathological. It exhibits large variation 
in the size and shape of its grains, as well as the number of neighbors each grain has. 
This data set allows us to explore the effectiveness of our approach in an extreme 
setting. In particular, it provides an example of a setting in which the standard 
choice of initial conditions is far from optimal. Finally, we consider an empirical 
data set: tomographic data obtained from a sample of Al-5 wt% Cu. We demon¬ 
strate that our method is able to produce an excellent approximation to this data 
set using a Laguerre tessellation. 


5.1. Artificial data 

The artificial data sets are produced using two distinct stochastic models. The 
polycrystalline model (PCM), developed in |8], describes a polycrystalline material 
undergoing grain coarsening. The second data set is produced using a randomly 
marked Poisson process; see m- The resulting tessellation is known as a Poisson- 
Laguerre tessellation (PLT); see |43j . In the first case, there is a strong correlation 
between the relative positions of the generating points and their weights. In the 
second example, the weights are independent of the positions of the generating 
points. As a result, the second data set is much less regular than the first data set. 
Figure [3] shows 2D cross sections of the data sets, together with their generating 
points. 

5.1.1. Model descriptions 

In the PCM model, a random Laguerre tessellation is produced using a set of spheres 
that may overlap slightly. The centers of the spheres give the locations of the gen¬ 
erating points and the radii define the weights. The spheres have hard-cores which 
may not overlap. This was shown to be a suitable model for polycrystalline ma¬ 
terials in J8] . Both the density of the sphere packing and the radii of the spheres 
influence the sizes and shapes of the cells in the resulting tessellation. In particular, 
when highly dense packings are used, very “spherical” cells are generated with a 
narrow coordination number distribution. We generated a sample from this model 
in a bounded window, W = [0, 700] 3 using parameters that were fitted to a sample 
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(a) Data set “artificial-PCM”. 


(b) Data set “artificial-PLT”. 


Figure 3. Partial cross-sections of artificial data sets. Spheres used as generators are drawn in gray (dots 
for sphere centers, gray circles for sphere surfaces). Radii are very large in the PLT case, therefore many 
surface voxels of spheres from other slices are visible. 


of Al-5 wt% Cu annealed for 200 minutes. The fitting procedure is described in (SJ. 
Note that the parameters we use were obtained by fitting the model to the empiri¬ 
cal data described in Section A2 Having produced a sample, we removed cells that 
would correspond to a very small number of voxels after discretization (i.e., cells 
having a smaller volume than a ball with radius of 4 voxels). The final result is a 
tessellation consisting of roughly 2500 grains (approximately the same number as 
in the experimental data). 

The PLT model is generated by simulating a homogeneous Poisson process with 
some intensity A in a bounded window W. The points are independently marked 
by gamma-distributed random variables with some shape parameter a and rate 
parameter f3. Note that, unlike the PCM model, it is possible that some generating 
points will not produce cells in the PLT model. We generated a sample from this 
model in W = [0, 700] 3 with intensity A = 7.5 x 1CD 6 (so that the expected number 
of points is roughly 2500), where we used parameters a = 50 and /3 = 1.125 for the 
mark distribution. 


5.1.2. Laguerre approximation 

The realizations of the PCM and PLT models were both discretized on a voxel grid, 
resulting in two 700 x 700 x 700 images. These images are of the form described 
in Section nm The CE method was then used to solve the LAP. The results are 
illustrated in Figure [4j where 2D cross-sections of the original tessellations are shown 
with the approximating tessellations superimposed. 

Using multi-threading on a standard quad-core processor (Intel Core i5-3570K), 
the computing time was roughly 3 hours for the PCM data and 4 hours for the PLT 
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data. For smaller test sets, with 500 grains each, the time required was roughly 
20 to 30 minutes. The memory requirements when fitting the full data sets were 
minimal (especially in comparison to approaches where it is necessary to store the 
complete voxelized data). Basically, it is necessary to store the test points and all 
the generators from one iteration of the CE method. The cost of storing additional 
variables, such as the means and standard deviations of the densities describing the 
generators and the cost values of the test points are negligible. For the PCM data, 
the number of test points used was approximately 140 000. Approximately 160 000 
points were used for the PLT data. Together with N ■ M ~ 2500 • 4000 generators in 
one iteration, this sums up to less than 200 MB of RAM when using 32-bit floating 
point coordinates. 



(a) PCM data. (b) PLT data. 


Figure 4. Partial cross-sections of artificial data: original cell boundaries (gray) superimposed with cell 
boundaries of our Laguerre approximation (black), with centers and radii of detected generators drawn in 
small gray dots and thin gray circles. 


Figure [5] illustrates the convergence behavior of the CE algorithm. When approx¬ 
imating the PCM data, roughly 650 iterations of the algorithm were required. In 
the PLT case, about 900 iterations were required. It is not surprising that the CE 


algorithm took longer to terminate in the PLT case. As will be seen in Section 5.1.3 
the initial conditions in this setting are much further from the optimal solution. 


5.1.3. Results 

In order to evaluate the quality of our approximations, we compare them to approx¬ 
imations obtained using the heuristic approach presented in [6] and the orthogonal 
regression method proposed in [08]. Note that, although the orthogonal regression 
approach is able to achieve quite good approximations of the cells, it does not result 
in a parametrized tessellation. Table [l] shows an evaluation of the approximations 
with respect to a volume-based discrepancy measure: the number of voxels that are 
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iteration in CE method 

(a) PCM data: cost values in CE iterations. 


iteration in CE method 

(b) PCM data: maximum of standard deviations 
in CE iterations. 
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200 400 600 


800 



iteration in CE method 

(c) PLT data: cost values in CE iterations. 


iteration in CE method 

(d) PLT data: maximum of standard deviations 
in CE iterations. 


Figure 5. CE method for artificial data: convergence of the cost values T> (Eq. §) and standard deviations. 


correctly labeled. Our method results in an almost perfect approximation of the 
PCM data. The heuristic approach considered in [6] also works very well. However, 
it is not able to reproduce the neighbor structure as successfully as our approach, 
as evidenced by rows 2 through 5 of Table [lj Our method is slightly less success¬ 
ful at approximating the PLT data but, in this case, considerably outperforms the 
heuristics of [6|. 

The parameters used by the heuristic approach of |6] are very close to the initial 
conditions of the CE algorithm. Thus, from the difference in performance, it seems 
that the initial conditions are quite far from the optimal parameters. In order to 
find better configurations, the CE algorithm has to be able to escape local minima 
around these initial conditions. This is made clear by the fact that increasing the 
number of times variance injection is used (as well as the number of iterations), we 
are able to further improve the results for the CE algorithm. Namely, by manually 
increasing the number of variance injections to 20 (instead of the 8 used in the 
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Table 1. Evaluation of the artificial data approximations: H denotes the heuristic approach 
of [6]; CE denotes the CE method considered in the present paper; OR denotes orthogonal 
regression proposed in 511- 



H 

PCM 

CE 

OR 

H 

PLT 

CE 

OR 

correctly labeled voxels [%] 

97.3 

99.1 

99.8 

87.8 

96.2 

99.8 

cells with all neighbors correct [%] 

75.6 

88.5 

78.9 

18.6 

45.1 

61.3 

cells with < 1 incorrect neighbors [%] 

96.3 

99.0 

97.3 

46.7 

77.3 

90.4 

cells with < 2 incorrect neighbors [%] 

99.6 

99.8 

99.7 

72.2 

92.1 

98.2 

mean number of erroneous neighbors/grain 

0.28 

0.13 

0.24 

1.83 

0.90 

0.51 

average displacement of centroids [voxels] 

0.56 

0.24 

0.09 

2.56 

1.33 

0.06 


standard approach), we obtained an approximation which correctly labeled 97.3% 
of the voxels (instead of 96.2%). This significant improvement indicates that lo¬ 
cal minima are present in which the CE algorithm is becoming trapped (because 
variance injection works by helping the algorithm escape local minima). This, in 
turn, implies that there are local minima near to the initial conditions of the CE 
algorithm that are sufficiently deep to trap it. A direct implication of these results 
is that algorithms that converge to nearby local minima (such as gradient-descent) 
may not lead to good approximations when using a standard choice of initial con¬ 
ditions. Furthermore, it is not clear that it is always straightforward to find good 
initial conditions. For example, the initial generating points are usually chosen to lie 
inside their cells. But, using the methods presented in m and the exact description 
of the PLT tessellation, it was not possible to find a solution to the LIP subject to 
the restriction that generating points were contained in their cells. If a method such 
as gradient-descent starts with all of the points lying inside their cells, it is unclear 
that it will be able to obtain a solution where some generating points lie far outside 
their cells. 

Figure [6] shows scatter plots of the original cell sizes vs. the cell sizes in the 
CE approximations. The quality of the CE approximation of the PCM data is 
immediately apparent. It is also clear that the difficulties in fitting the PLT data 
lie primarily in the small cells. 


5.2. Experimental data 


The experimental data we consider was used in [8j to develop the PCM model 


discussed in Section 5.1 An Al-5 wt% Cu sample (cylindrical with 8.5 mm length 


and 4 mm diameter) was heated to the semisolid state at a temperature of 592 °C, 
at which point coarsening processes were measured in situ using synchrotron X-ray 
tomographic imaging. For the present paper, we consider only data obtained after 
an annealing time of 200 minutes (which corresponds to the first annealing step 
in ED- It consists of approximately 2500 grains, which is quite a large data set. A 
2 D cross-section of the data is shown in Figure [ 7 J The experimental data set was 
segmented using the watershed transformation, resulting in a labeled image of grain 
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cell size in original data 

(a) PCM data. 


cell size in original data 

(b) PLT data. 


Figure 6. Scatter plot of cell sizes given as radii of volume-equivalent spheres obtained from the original 
data and cells extracted by the CE method. 


regions, which are separated by a watershed layer of one-voxel thickness. Further 
details regarding sample preparation, imaging and segmentation can be found in 

0 . 



(a) Tomographic grayscale data. (b) Grain boundaries obtained by watershed 

segmentation. 

Figure 7. Partial cross-section of experimental data. 


5.2.1. Laguerre approximation 

The CE method was applied to solve the LAP for the experimental data set. A 
cross-section of the resulting Laguerre approximation is given in Figure [8} 

Fitting the approximation took approximately 70 minutes on the same Intel Core 
i5-3570K quad-core processor used to fit the artificial data. The number of test 
points extracted from the image data was roughly 150 000, which is almost the 
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Figure 8. Partial cross-section of experimental data: original boundaries of grain regions (gray) super¬ 
imposed with cell boundaries of our Laguerre approximation (black), with centers and radii of detected 
generators drawn in small gray dots and thin gray circles. 


same as for the artificial data sets. The convergence behavior of the CE method 
is illustrated in Figure [9j Using the same parameters, only 283 iterations were 
necessary, which is substantially fewer than the number required for the artificial 
data. This is mainly due to the fact that only one variance injection was required. 




iteration in CE method iteration in CE method 

(a) Cost values for experimental data in CE (b) Maximum of standard deviations for 
iterations. experimental data in CE iterations. 

Figure 9. CE method for experimental data: convergence of the cost values T> (Eq. and standard 
deviations. 


5.2.2. Results 

As with the artificial data, we compare our results with the heuristic approach pro¬ 
posed in [6| and the orthogonal regression approach introduced in |48]. As mentioned 
above, the orthogonal regression approach reconstructs the individual cells quite well 
but does not result in a tessellation. The results are summarized in Table [2j The CE 
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method correctly assigns 91.4% of the voxels. In contrast, the heuristic approach 
correctly assigns only 82.9% of the voxels. The orthogonal regression yields 96.3%. 


Table 2. Evaluation of the experimental data approximations: H de¬ 
notes the heuristic approach of [6]: CE denotes the CE method consid¬ 
ered in the present paper; OR denotes orthogonal regression proposed 
in [48]. 



H 

CE 

OR 

correctly labeled voxels [%] 

82.9 

91.4 

96.3 

cells with all neighbors correct [%] 

20.1 

38.4 

29.3 

cells with < 1 incorrect neighbors [%] 

48.4 

73.1 

50.2 

cells with < 2 incorrect neighbors [%] 

65.6 

90.6 

62.6 

mean number of erroneous neighbors/grain 

2.29 

1.04 

2.37 

average displacement of centroids [voxels] 

3.40 

1.78 

0.23 


Note that the CE method significantly outperforms both the heuristic method and 
orthogonal regression when describing the neighborhood structure of the grains, al¬ 
though it does not seem possible for a tessellation containing only convex cells to 
accurately capture the full neighborhood structure. This is at least partially due 
to segmentation issues and contacts in the image data with very small areas. In 
contrast, the orthogonal regression approach performs surprisingly badly. It seems 
that the geometric properties of normal Laguerre tessellations (e.g., coinciding faces, 
edges and vertices) favor realistic reconstructions. We think this confirms that the 
Laguerre tessellation is a good choice for representing polycrystalline microstruc- 
tures. It is possible, of course, that the quality of the fit could be improved using 
non-convex cells. For example, in [Mj, the linear programming method was used 
to obtain a generalized power diagram approximation of similar data (but with far 
fewer grains and voxels) resulting in a fit that correctly assigned 93.8 % of the voxels 
with non-convex cells. 



cell size in original data cell size in original data 

(a) Scatter plot of cell sizes obtained from (b) Distances between centroids of original and 
original and extracted cells. approximated cells, in dependence on cell sizes. 


Figure 10. Experimental data: comparison of original data and cells extracted by the CE method. Cell 
sizes are given as the radii of their volume-equivalent spheres. 
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Figure [To|a) shows a scatter plot of the volumes of the original grains against 
the volumes of the Laguerre cells (with the volumes expressed as radii of volume- 
equivalent spheres). The overall fit of the cell volumes is excellent. Figure 10 b) 
shows that the locations of the grains are also quite accurate. Note that the main 
issue with fitting seems to be small cells. This is not such a problem, however, as 
small cells (with equivalent radii of up to 10 voxels) are subject to image segmen¬ 
tation error and are, thus, not too reliable in the original data. 


6. Conclusions and outlook 

In this paper, we considered the problem of approximating tomographic data by 
a Laguerre tessellation. We expressed this problem as an optimization problem: 
the generating points of the approximating tessellation need to be chosen in order 
to minimize the discrepancy between the tessellation and the tomographic data. 
We considered an interface-based discrepancy measure, instead of the volume-based 
discrepancies more commonly considered in the literature. This allowed us to use the 
CE method, a stochastic optimization method that is able to escape local minima, 
to fit the approximating tessellation. We then carried out numerical experiments 
on both artificially generated and experimentally obtained tomographic data that 
demonstrated the broad effectiveness of our approach. 

Our method is robust and easy to implement. Thus it can be applied to fit La¬ 
guerre tessellations to almost any material with a granular or cellular structure. An 
obvious next step is to extend our approach to tessellations that include non-convex 
cells. In particular, we believe the approach and philosophy outlined in this paper 
can be extended to fit generalized power diagrams (see m ED) — generalizations 
of Laguerre tessellations that can include non-convex cells. This will be the subject 
of a forthcoming research paper. 
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Supplemental material 

A software package including Java code and data sets can be downloaded from 

https://github.com/stochastics-ulm-university/laguerre-approximation. 
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